perm filename APPEND.PUB[NSF,MUS]1 blob sn#096512 filedate 1974-04-10 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00009 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	.SELECT B
C00014 00003	.GROUP SKIP 1
C00017 00004	.NEXT PAGE
C00023 00005	.GROUP SKIP 2
C00028 00006	.NEXT PAGE
C00031 00007	.NEXT PAGE
C00035 00008	.GROUP SKIP 1
C00038 00009	.SELECT A
C00045 ENDMK
C⊗;
.SELECT B
.BEGIN CENTER
III. APPENDICES

.SELECT A
APPENDIX A: THE HETERODYNE FILTER
.END
.SELECT 1
.BEGIN ADJUST FILL
As was mentioned before,
in additive synthesis we physically model a complex sound waveform as
a  sum of sinusoids  with slowly time-varying  amplitudes and phases.
The process of synthesis involves specifying the amplitude  and phase
for each  component
sinusoid as  it varies with time throughout the duration of the tone.
These sinusoids  are  added
together to  produce the  complex waveform.  Equation (1) is
duplicated here as equation (A1)  and summarizes
this  formulation.
The  sinusoids  in equation  (A1)  represent  the
components of the complex tone.
.END
.GROUP
.SELECT 3


          M
(A1) F%8α%3 = %6S%3 A%8n%3 sin(%4w%8n%3αh+%4q%8n%3)
         n=1

Notation:    h is the time between consecutive samples
	     α is the sample number
	     F%8α%3 is the sampled, digitized waveform at time αh
	     A%8n%3 is the amplitude of the nth partial tone
	         and is assumed to be slowly varying with time
	     %4q%8n%3 is the phase of the nth partial tone
	         and is assumed to be slowly varying with time
	     %4w%8n%3 is the radian frequency of the nth partial tone
		

.SELECT 1
.APART
.BEGIN FILL ADJUST
To be more explicit, we should point out again that A%8n%1 and %4q%8n%1
are functions of time. We will indicate this in the following text by
appending the subscript %3α%1 to each of these functions:
.GROUP
.SELECT 3

            A%8n%3 = A%8nα%3  and  %4q%8n%3 = %4q%8nα%3

.SELECT 1
.APART
This brings out the time dependence more clearly.

It should be pointed out that this is not a Fourier series representation
despite its outward similarity. The Fourier series approximates a periodic
function with a sum of orthogonal sinusoids with fixed frequencies and
phases and (possibly) exponential amplitudes. The model described in equation
(A1) is not necessarily periodic, has time varying phase, and has time varying
amplitudes that are not necessarily exponential. The individual componants
are certainly not orthogonal in the general case.

We have avoided use of continuous analysis and have defined all our
functions in the discrete domain. This is done because the ultimate realization
of these processes is on a digital computer. Inaccuracies can result
from doing the mathematics in the continuous case and assuming it
can be converted to the discrete domain merely by sampling. The conversion
to the discrete domain must be done with some care.
Since discrete mathematics and the science of digital signal processing
have become so well developed, there is little reason to define our
processes in the continuous domain, only to subsequently realize them in the
discrete domain.

We have borrowed the notation of numerical analysis for functions sampled
at equally spaced intervals by placing the sample number as a subscript.
This is done partly for notational convenience and partly because in the
computer, sampled functions are represented by matrices.

We seek a method of determining the
functions A%8nα%1 and %4q%8nα%1 of a tone from a musical instrument,
we can then synthesize an approximation to the waveform F%8α%1 from those functions
by use of equation (A1).

We write equation (A1) using sinusoids at time-varying phase angles rather
than as sine and cosine quadrature componants because it is more efficient
to synthesize the waveform in the form shown in equation (A1).
.END
.GROUP SKIP 2

.BEGIN FILL ADJUST
Let us turn to the problem of determining the functions
A%8nα%1 and %4q%8nα%1 of a music   instrument tone. To aid the
analysis, we must assume the frequencies of the partial tones,
%4w%8n%1, are nearly harmonically related. That is, there is some
frequency, %4w%1, such that %4w%8n%1 is approximately n%4w%1. We shall
call this frequency %4w%1 the fundamental frequency of the tone.

Basically, the method is as follows:
First, compute the following two summations at each point in time %3α.
.END
.SELECT 3

         α+N-1
(A2) %9a%8n%8α%3 = %6S%3 F%8i%3sin(n%4w%80%3ih+%4f%80%3)
          i=α

         α+N-1
(A3) %9b%8n%8α%3 = %6S%3 F%8i%3cos(n%4w%80%3ih+%4f%80%3)
          i=α

%1
.BEGIN FILL ADJUST
The initial phase angle, %4f%80%1 is included for generality.
It will be shown below that the method is independent of the
initial phase angle.
From these, we calculate two more sequences:
.END

%3(A4)    A%8nα%3 = (%9a%8nα%22%3+%9b%8nα%22%3)%21/2%3

(A5)    %4q%8nα%3 = atan(%9a%8nα%3/%9b%8nα%3)

.SELECT 1
.BEGIN FILL ADJUST
The summations are taken to be over one period of a sinusoid
of frequency %4w%80%1, that is, Nh%4w%80%1 = 2π. This places
somewhat of a restriction on the frequency of analysis, %4w%80%1,
because in the discrete domain, the period, N, is restricted
to integral values. This has not proved a problem in our experience.

If the partial tones are nearly harmonically related, if the
amplitude and phase functions of the tone vary slowly with time,
and if %4w%80%1 is nearly
equal to the fundamental frequency of the tone, then A%8nα%1 and
%4q%8nα%1, as computed by equations (A2) through (A5),
will indeed be approximations to the actual amplitudes and
phases of the partials of the tone under analysis.

To review, we do the computations indicated in
equations (A2), (A3), (A4), and (A5) for each of the partials of a
tone, over the entire time interval spanned by the tone. The output
A%8nα%1 and %4q%8nα%1 may then be used in equation (A1) to synthesize a
tone based on the analysis.
.END
.GROUP SKIP 2
.BEGIN ADJUST FILL
Now let us compute the response of the heterodyne
filter as defined by equations (A2), (A3), (A4), and (A5) to a
sinusoid of constant amplitude and phase. We do this by substituting
for F%8i%1 in equations (A2) and (A3) the function sin(%4w%1ih).
We may compute the summations without error by use of the summation
calculus (Hamming). Using the fact that Nh%4w%80%1 = 2π, we may
calculate A%8nα%1 and %4q%8nα%1 explicitly.
.END
.GROUP
.SELECT 3

            1                     1                 1
(A6)  A%8nα%3 = --- sin%22%3(%4w%3Nh/2)%9{%3----------------- + -----------------
           4N%22%3              sin%22%3[(%4w%3+n%4w%80%3)h/2]   sin%22%3[(%4w%3-n%4w%80%3)h/2]



                            2cos[n%4w%80%3h-2%4f%80%3]
                  + --------------------------------%9}%3
                    sin[(%4w%3+n%4w%80%3)h/2] sin[(%4w%3-n%4w%80%3)h/2]


.APART
.SELECT 1
.BEGIN FILL ADJUST
The expression for %4q%8nα%1 is a bit long and is thus not included here.
As we consider the limit as %4w%1 approaches n%4w%80%1, we find great
simplification of the results. Let us define %4Dw%1 as (%4w%1-n%4w%80%1).
.END
.GROUP
.SELECT 3

                   1
(A7)   lim   A%8nα%3 = ---%9{%30 + N%22%3 + 0%9}%3 = 1/4
       %4w%3→%4w%80%3       4N%22%3


                    sin%9{%32%4w%80%3h[(N-1)/2+α]%9}%3 + N sin%9{%4Dw%3h[(N-1)/2+α]%9}%3
(A8)   lim %9a%8nα%3/%9b%8nα%3 = ----------------------------------------------
       %4w%3→%4w%80%3         cos%9{%32%4w%80%3h[(N-1)/2+α]%9}%3 + N cos%9{%4Dw%3h[(N-1)/2+α]%9}%3


.APART
.NEXT PAGE
.SELECT 1
.GROUP
If N>>1 then (A8) reduces greatly.
.SELECT 3


(A9)   lim %9a%8nα%3/%9b%8nα%3 = tan%9{%4Dw%3h[(N-1)/2+α]%9}%3
       %4w%3→%4w%80%3


.APART
.GROUP SKIP 1
.SELECT 1
.BEGIN FILL ADJUST
Thus we see that in the limit, A%8nα%1 approaches one quarter of the
amplitude of the input sinusoid and %4q%8nα%1 is a term related to
the difference in frequency of the input sinusoid with the analysis
frequency. With instruments with partials whose frequencies deviate
from the ideal, this provides a dynamic estimation of those frequencies.

Notice that neither A%8nα%1 nor %4q%8nα%1 are functions of %4f%80%1, the
initial phase angle. Under the assumptions stated,
this process is independent of initial phase.

Notice also that A%8nα%1 is no longer a function of %3α%1, the time parameter,
but %4q%8nα%1 is. Fortunately, %4q%8nα%1 is a linear function of %3α%1
whose slope is simply %4Dw%1h.

Figure A1 shows a plot of A%8nα%1 as indicated in equation (A1) for a
range of values of %4w%1. In this case, %4w%80%1 is 2π(125 Hz) and
n is 4. We see that there is a zero of transmission at all integral
multiples of %4w%80%1 except the n%2th%1 multiple.

This technique is useful as long as the amplitudes and phases
of the partials of the input waveform change  slowly with time. If
the frequencies of the partials deviate from integral multiples of the fundamental
by too great an amount, further error may be introduced.

One must remember that the heterodyne filter as defined here is
a %5nonlinear%1 filter. Although we derived the response to a sinusoid,
superposition does not hold in such a filter. The calculation of
%9a%8nα%1 and %9b%8nα%1, however, is quite linear and superposition again
applies. The general statements made above still hold, but only because of the
special nature of the input signal.
.END
.GROUP SKIP 2
.NEXT PAGE
.BEGIN CENTER
.SELECT A
APPENDIX B: MULTIDIMENSIONAL SCALING TECHNIQUES
.END
.GROUP SKIP 1
.SELECT C
SPATIAL MODEL FOR PERCEPTUAL RELATIONSHIPS
.SELECT 1
.BEGIN FILL ADJUST
It has been found most useful to employ  a spatial model to represent
the judged  relationships between sets of stimuli,   such as auditory
signals.   Computer-based  multidimensional scaling  algorithms  have
been developed  for the reduction  of the very complex  data obtained
from  the subjective  evaluations of  perceived relationships between
all members in a  set of stimuli. The  result is displayed in  a form
which  is  much  more  easily  comprehended and  interpreted  by  the
investigator, that  of  a  geometric configuration  of  points  which
represent the  individual stimuli.   The structure of  the subjective
evaluations   of  the  set  of   stimuli  is  then   mapped  into  an
n-dimensional space,   where  the distances  between  the points  are
determined by some measure of  the psychological distance between all
pairs of stimuli.

The  psychological measures which can be  mapped into a spatial model
in  terms  of  distance  include  the  confusability  or  the  judged
similarity  of  stimuli.   One  experimental  procedure involves  the
identification  of  individual   signals,  perhaps  presented   under
different conditions, and the result is a  square confusion matrix of
stimulus by response.  A transform of this matrix yields the relative
psychological distances, directly related  to confusability,  of  all
points in  the data  matrix. Another experimental  procedure involves
the  rating of relative similarity for a  pair of stimuli in the set.
The similarity ratings can be placed in another  square matrix, where
each entry  point is for one  particular of all  the possible ordered
pairs of  stimuli.    The  psychological distance  in  this  case  is
inversely  related to  the  similarity of  stimuli,  and is  directly
related to their dissimilarity.

The common  characteristic of the scaling programs  we find useful is
their  generation  of  an  empirically-based  representation  of  the
relationships between  the stimuli,   rather than  some theoretically
imposed, a priori representation. We proceed from the perceptual data
and  will compare  the  representation  of  this data  to  the  known
physical attributes of the stimuli.  The uncovering of psychophysical
relationships, that is, making  correlations between the  subjective,
psychological judgments of the stimuli and their physical properties,
is  essentially a  matter of interpreting  the representation  of the
perceptual data in terms of the known physical data for  the stimuli.
Note  the  significant  difference  of this  from  most  experimental
approaches,  which begin with some a priori  notions of the nature of
the results,  and the  experiment  is constructed  and data  analyzed
around these notions. The multidimensional scaling approach is useful
when the stimuli  are inherently so  complex that we  have no a
priori notions and are  willing to rely totally on  the empirical
method of analysis.   We are concerned with both the dimensionality and
the general  properties of  the space.    Correlations with  physical
parameters are  sought.  Various  programs exist which are  useful in
assessing  the correspondence of  data structures,  so  that it would
prove fruitful to formally represent the physical data.
.END
.GROUP SKIP 2
.SELECT C
MULTIDIMENSIONAL SCALING ALGORITHMS
.SELECT 1
.BEGIN FILL ADJUST
We  will  briefly  describe the  two  programs  for  multidimensional
scaling found most  useful in our research, MDSCAL and INDSCAL. These
programs both attempt to represent input data matrices in the form of
a  configuration  of points  located  in  an n-dimensional  geometric
space, where n, the number of dimensions, is specifiable by the user.
The  points correspond to the stimuli,  whose psychological distances
are given in the input matrices.  The coordinates  of the  points are
obtained by an iterative computational algorithm which optimizes  the
correspondence of interpoint  distances in the spatial representation
to the measured psychological distances between the stimuli. 

  
MDSCAL peforms a  non-metric multidimensional  scaling.  The  optimal
spatial  representation of  a subjective  response  matrix is  one in
which the rank order of  the values of psychological distance be  the
same  as  the   rank  order  of  the  interpoint   distances  in  the
n-dimensional  geometric  configuration.  A  monotonic  function maps
psychological values  into distances  in  the spatial  representation
(for a discussion of the theory and procedures of this algorithm  see
Shepard, 1962a, 1962b; Kruskal  1964a, 1964b).  The use  of MDSCAL in
conjunction  with   other  programs  which  deal  with  psychological
distance matrices can be particularly informative.  One such program,
HICLUS (Johnson, 1967) produces a tree-structure which represents the
hierarchical clustering  relationships of  stimuli in  the matrix  as
inferred from their  psychological distances (the use  of HICLUS with
MDSCAL  for the analysis of confusion  matrices for speech signals is
demonstrated by Shepard, 1972).
.NEXT PAGE
INDSCAL  is a metric  multidimensional scaling  program,   which  was
developed to  utilize the individual differences  in sets of response
matrices  for  analysis.     It  generates  a   single  n-dimensional
representation  for the complete  set of  matrices.  It  analyzes the
variations  in the  set  of  individual  data  matrices  to  uniquely
determine a  rotation for the  axes in the  space.  Also  produced by
this  analysis is  a representation  of weightings which  account for
individual response  variations  along  the spatial  dimensions.  The
weightings  are mapped  into an  n-dimensional  spatial configuration
which can  be used  to assess  the relationships  between  individual
subjects  or  experimental  conditions  (see  Carroll,  1970,  for  a
complete discussion  of the theory and  operations of this technique;
and Carroll and Wish, 1973,  for an application of the method for  the
representation of confusion matrices for speech signals).
.END
.NEXT PAGE
.SELECT A
APPENDIX C: SPECTRAL SHAPING FILTERS
%1
.BEGIN FILL ADJUST
We use several basic filters. We shall describe only
the ones that are different from the more common low-pass
and high-pass networks, as these are documented extensively
in the literature.
.END
.SKIP 1
.SELECT C
SIMPLE RESONATORS AND ANTI-RESONATORS
.SELECT 1
.BEGIN FILL ADJUST

Let us begin with the digital resonator. This is most simply
expressed as follows:
.END
.SELECT 3

                     1-qZ%2-1%3
(C1)    T(Z) = ----------------------
               1-2r cos(%4w%80%3h)Z%2-1%3+r%22%3Z%2-2%1


With recurrence relation%3

(C2)	Y%8n%3 = X%8n%3-qX%8n-1%3+2r cos(%4w%80%3h)Y%8n-1%3-r%22%3Y%8n-2%1

.BEGIN FILL ADJUST
This has a resonant frequency of %4w%80%1. r determines the Q of the
resonator. q is a zero of transmission. This resonator is used
extensively by Gold and Rader (1969), and thus will not be examined
further here.

Another filter that is useful is the anti-resonance,
or notch filter. In the continuous case, the transfer function
is given by:
.END
%3
               S%22%3+2%4d%3S+%4d%22%3+%4w%80%22%3
(C3)    H(S) = ------------
               S%22%3+2%4s%3S+%4s%22%3+%4w%80%22%1

.BEGIN FILL ADJUST
To prevent warping of the frequency axis, the Boxer-Thaler (1956)
transformation is used :
.END
.GROUP
%3

                 C%81%3+C%82%3Z%2-1%3+C%83%3Z%2-2%3
(C4)    T(Z) = --------------
                 C%84%3+C%85%3Z%2-1%3+C%86%3Z%2-2%1


.APART
.NEXT PAGE
.GROUP
With its related recurrence formula:

%3(C5)	Y%8%3 = (C%81%3X%8n%3+C%82%3X%8n-1%3+C%83%3X%8n-2%3-C%85%3Y%8n-1%3-C%86%3Y%8n-2%3)/C%84%3

    where  C%81%3 = h%22%3(%4d%22%3+%4w%80%22%3)+12h%4d%3+12
           C%82%3 = 10h%22%3(%4d%22%3+%4w%80%22%3)-24
           C%83%3 = h%22%3(%4d%22%3+%4w%80%22%3)-12h%4d%3+12
           C%84%3 = h%22%3(%4s%22%3+%4w%80%22%3)+12h%4s%3+12
           C%85%3 = 10h%22%3(%4s%22%3+%4w%80%22%3)-24
           C%86%3 = h%22%3(%4s%22%3+%4w%80%22%3)-12h%4s%3+12%1

.APART
.BEGIN FILL ADJUST
This gives a magnitude-frequency response that is unity everywhere
except near %4w%80%1. The strength at %4w%80%1 is determined by
the values of %4s%1 and %4d%1. If both %4s%1 and %4d%1 are small
compared to %4w%80%1, then the magnitude of H(S) at S=j%4w%80%1
will be approximately %4d%1/%4s%1. If %4d%1 is zero, the magnitude
goes to identically zero at %4w%80%1. If %4d%1=%4s%1, then we have
an all-pass network like the one used in appendix B as the
second order reverberator.

The magnitude-frequency response of the filter in equation (C4) is
plotted in figure C1 for five different values of %4d%1/%4s%1. The
value of %4w%80%1 was 2π(500 Hz), and %4s%1 was 2π(50 Hz). The
ratios %4d%1/%4s%1 were 4, 2, 1, .5, and .25.
.END
.GROUP SKIP 1
.SELECT C
THE COMB FILTERS
.SELECT 1
.BEGIN FILL ADJUST
The comb filter comes in four forms, two of which are all-zero filters
and two of which are all-pole filters. The two all-zero filters
have the following recurrence relations:
.END
.GROUP

%3(C6)    Y%8n%3 = X%8n%3+X%8n-m%3

(C7)    Y%8n%3 = X%8n%3-X%8n-m%1

.APART
.BEGIN FILL ADJUST
These have magnitude frequency responses as follows:
.END

%3(C8)    |T(%9e%3↑j%4↑w%3↑h)| = {[1+cos(m%4w%3h)]%22%3+sin%22%3(m%4w%3h)}%21/2%3

(C9)    |T(%9e%3↑j%4↑w%3↑h)| = {[1-cos(m%4w%3h)]%22%3+sin%22%3(m%4w%3h)}%21/2%1

.NEXT PAGE
.BEGIN FILL ADJUST
It is clear that equation (C8) is zero at m%4w%1h = (2n+1)π and equation
(C9) is zero at m%4w%1h = 2nπ. Thus, applying the filter specified by
equation (C7) to a perfectly periodic waveform of period mh/n
of arbitrary harmonic content will exactly annihilate the waveform. The
filter of equation (C6) is somewhat more subtle. It will annihilate, for instance,
the odd harmonics of a waveform of period 2mh.
.END
.GROUP SKIP 1
.BEGIN FILL ADJUST
By placing a constant in the recurrence relations, we may move the
zeros off the j%4w%1 axis:
.END
.GROUP

%3(C10)    Y%8n%3 = X%8n%3+gX%8n-m%3

(C11)    Y%8n%3 = X%8n%3-gX%8n-m%1

.APART
.BEGIN FILL ADJUST
This causes the frequency response to approach zero, but never become
identically zero unless g is exactly unity.

We may invert the spectrum of these filters by placing the delay
in the feedback path, making recursive filters of these two:
.END

%3(C12)   Y%8n%3 = X%8n%3-gY%8n-m%3

(C13)   Y%8n%3 = X%8n%3+gY%8n-m%1

.BEGIN FILL ADJUST
The fact that the signs are reversed in the two equations is no
accident. Equation (C12) does indeed correspond to equation (C10),
and likewise equation (C13) to (C11). These two filters are like the
previous ones except that they have resonances where the others
had anti-resonances.

The magnitude-frequency responses for the filters in equations (C7)
and C(11) were plotted in figures C2 and C3 respectively
for 1/(mh) = 400 Hz and for four different values
of g, which were .25, .5, .75, and 1. It might be noted that the
graphs are truncated at 20 db. Actually, for g=1, the response in figure
C2 goes to zero, which would require our plot to extend to -%3∞%1. Likewise,
the response in figure C3 goes to +%3∞%1.
We have somewhat arbitrarily placed a maximum excursion on the plot at
+ or - 20 db.
.END
.NEXT PAGE
.SELECT A
APPENDIX D: ALL-PASS UNIT REVERBERATORS
.SELECT 1
.GROUP SKIP 1
.BEGIN FILL ADJUST
A particularly useful unit of artificial reverberation is the all-pass network.
We use two different forms:
(1) a first order unit of which the impulse response is a pulse train with
exponentially decaying amplitude, and (2) a second order unit
of which the impulse response is a pulse train of which the amplitude is
a damped sinusoid.

A simple all-pass filter is the pole-zero pair, symmetrically
located about the j%4w%1 axis. The complex frequency response in the
continuous case is shown in equation (D1).
.END

%3
	       S-%4s%3
(D1)	H(S) = ----
	       S+%4s
%1
.BEGIN FILL ADJUST
where S is complex frequency and 1/%4s%1 is the decay time. This filter
has a pole at -%4s%1 and a zero at %4s%1.
We may convert this to a digital filter by use of the bilinear transform
(Gold and Rader).
.END

%3
                 (%4s%3h+2)Z%2-1%3 + (%4s%3h-2)
(D2)    T(Z) = - -----------------
                 (%4s%3h-2)Z%2-1%3 + (%4s%3h+2)
%1
.BEGIN FILL ADJUST
where Z%2-1%1 is the unit delay operator and h is the time between
consecutive samples.
Since this is an all-pass, the magnitude of T(Z) anywhere on the unit
circle is the same. This being the case, we may raise Z in equation (D2) to
any power without altering the magnitude of T(Z) on the unit circle.
If we raise Z to the power of -m, the frequency response will cycle m
times as we go around the unit circle.
The decay time becomes m/%4s%1. 
.END

%3
                 (%4s%3h+2)Z%2-m%3 + (%4s%3h-2)
(D3)    T(Z) = - -----------------
                 (%4s%3h-2)Z%2-m%3 + (%4s%3h+2)
%1

Which implies the recurrence relation:

%3
(D4)	Y%8n%3 = {(%4s%3h-2)X%8n%3+(%4s%3h+2)X%8n-m%3-(%4s%3h-2)Y%8n-m%3}/(%4s%3h+2)
%1
.BEGIN FILL ADJUST
This is essentially the unit reverberator used by Schroeder (1961),
except we have realized it in the canonical form here, thus saving one
multiplication over the form previously used. The frequency
response is identically constant around the unit circle.
The impulse response is a pulse train with exponentially decaying amplitude.
The ratio (%4s%1h-2)/(%4s%1h+2) is called the gain of the reverberator.
As the gain approaches unity, the impulse response of the reverberator
decays more and more slowly.
We call this the `first-order' unit reverberator, although technically
it is of order m.
A block diagram of this unit and a plot of a typical
impulse response is shown in Figures 13a and 13b.

There is one important generalization of the all-pass network.
If we begin with an all-pass filter which has complex conjugate
poles rather than a single real pole, we realize a reverberator that
differs significantly in character from the one in equation (D3).
We begin with the following filter, again shown in the continuous
case:
.END

%3
               (S%22%3-2%4s%3S+%4s%22%3+%4w%80%22%3)
(D5)    H(S) = ---------------
               (S%22%3+2%4s%3S+%4s%22%3+%4w%80%22%3)
%1
.BEGIN FILL ADJUST
where %4w%80%1 is the resonant frequency of the filter.
This is again an all-pass. Let us transform it via the bilinear
transform and substitute a unit delay of m as was done above.
.END

%3
               C%81%3Z%2-2m%3+C%82%3Z%2-m%3+C%83%3
(D6)    T(Z) = --------------
               C%83%3Z%2-2m%3+C%82%3Z%2-m%3+C%81%3
%1

Which, in turn, implies the recurrence relation:

%3
(D7)	Y%8n%3 = (C%83%3X%8n%3+C%82%3X%8n-m%3+C%81%3X%8n-2m%3-C%82%3Y%8n-m%3-C%83%3Y%8n-2m%3)/C%81%3

    where  C%81%3 = %4w%80%22%3h%22%3+%4s%22%3h%22%3+4%4s%3h+4
           C%82%3 = 2%4w%80%22%3h%22%3+2%4s%22%3h%22%3-8
           C%83%3 = %4w%80%22%3h%22%3+%4s%22%3h%22%3-4%4s%3h+4
%1
.BEGIN FILL ADJUST
This, then, is another all-pass unit reverberator of a different character,
where the impulse response is a pulse train of which the amplitude
is a damped sinusoid.  We call this the `second-order' unit reverberator.
It should be noted that with the substitution of Z%2-m%1
for Z%2-1%1, the frequency of the sinusoid for this unit
(equation (D6)) becomes %4w%80/m%1. The decay time likewise becomes m/%4s%1.

A block diagram of this unit and a plot of a typical impulse response is
shown in figures 14a and 14b.
.END